clear;

delta_exac = 0;
delta_predict = 0;
cnt = 0;
filename = "case9";

for i = 1:1000
    mpc = loadcase(filename);
    mpc = ext2int(mpc);
    mpc.bus(:,3) = 0;
    mpc.gen(:,2) = 0;
    
    gen_buses = find(mpc.bus(:, 2) == 2 | mpc.bus(:, 2) == 3);
    load_buses = find(mpc.bus(:, 2) == 1);
    %randomize the raw data
    gen_selected = gen_buses(randperm(length(gen_buses), round(0.3 * length(gen_buses))));
    gen_selected_idx = find(ismember( mpc.gen(:,1), gen_selected));
    load_selected = load_buses(randperm(length(load_buses), round(0.3 * length(load_buses))));
    
    a = normrnd(0, 0.5);
    mpc.bus(load_selected, 3) = (1+a) * mpc.bus(load_selected, 3);
    mpc.bus(load_selected, 4) = (1+a) * mpc.bus(load_selected, 4);
    b = normrnd(0, 0.3);
    mpc.gen(gen_selected_idx,6) = (1+b) * mpc.gen(gen_selected_idx,6);
    
    %predict whehter the system is stable
    Ybus = full(makeYbus(mpc));
    Ybus = imag(Ybus);
    
    BGG = Ybus(gen_buses, gen_buses);
    BLL = Ybus(load_buses, load_buses);
    BLG = Ybus(load_buses, gen_buses);
    [~, genidx] = ismember(gen_buses, mpc.gen(:,1));
    VG =  mpc.gen(genidx,6);
    V_L_star = -inv(BLL) * BLG * VG;
    Q_crit = 0.25 * diag(V_L_star) * BLL * diag(V_L_star);
    QL = mpc.bus(load_buses, 4)/mpc.baseMVA;
    temp = Q_crit\QL;
    Delta = norm(temp, "inf");
    if Delta <= 1
        suc = 1;
    else
        suc = 0;
    end
    
    %exact solution
    mpopt = mpoption('out.all', 0, 'verbose', 0);
    [result, success] = runpf(mpc, mpopt);
    
    if success == 0
        continue;
    else
        bus_voltages = result.bus(load_buses, 8);
    end
    
    %whether the prediction is right or not
    if success == suc
        cnt = cnt + 1;
        delta_predict = delta_predict + (1 - sqrt(1 - Delta))/2;
        delta_exac = delta_exac + norm(abs(bus_voltages - V_L_star)./V_L_star,"inf");
    else
        disp("It fails.")
    end
end
accuracy = cnt / 10;
delta_exac = delta_exac / cnt;
delta_predict = delta_predict / cnt;
condition_accuracy = abs(delta_exac - delta_predict)/delta_exac;


accuracy
delta_predict
delta_exac
condition_accuracy

